---
title: "RACIAL BIAS AND PUBLIC SUPPORT FOR US DRONE STRIKES = Second Survey Experiment Replicaiton Kit"
author: "Paul Lushenko"
date: "6/26/2024"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
library(plyr)
library(dplyr)
library(ggplot2)
library(readr)
library(ggpubr)
library(tidyverse)
library(tidyr)
library(stargazer)
library(corrplot)
library(Hmisc)
library(margins)
library(magrittr)
library(wordcloud)
library(tm)
library(grid)
library(coefplot)
library(gplots)
library(car)
library(rstatix)
library(emmeans)
library(lmtest)
library(sandwich)
library(effects)
library(Interact)
library(sciplot)
library(lsr)
library(qwraps2)
library(gtsummary)
library(kableExtra)
library(caret)
library(ggthemes)
library(ggprism)
library(patchwork)
library(magrittr)
library(rstatix)
library(sjmisc)
library(scales)  
library(modelsummary)
```

```{r read in data, include=FALSE}
race <- read.csv("",
          stringsAsFactors = FALSE,                                                           
          na.string = ".") %>% filter(Finished == 1 & participate == 2) #for main analysis
```

```{r format data, include=FALSE} 
race <- race[-(1:2),]
names <- as.data.frame(names(race))
race <- race %>% 
  select(-`StartDate`,
         -`EndDate`,
         -`Status`,
         -`IPAddress`,
         -`Progress`,
         -`Finished`,
         -`RecordedDate`,
         -`ResponseId`,
         -`RecipientLastName`,
         -`RecipientFirstName`,
         -`RecipientEmail`,
         -`ExternalReference`,
         -`LocationLatitude`,
         -`LocationLongitude`,
         -`DistributionChannel`,
         -`UserLanguage`,
         -`Q_RecaptchaScore`,
         -`participate`,
         -`Q6A_First.Click`,
         -`Q6A_Last.Click`,
         -`Q6A_Click.Count`,
         -`Q7A_First.Click`,
         -`Q7A_Last.Click`,
         -`Q7A_Click.Count`,
         -`Q8A_First.Click`,
         -`Q8A_Last.Click`,
         -`Q8A_Click.Count`,
         -`rid`,
         -`age`,
         -`gender`,
         -`hhi`,
         -`ethnicity`,
         -`hispanic`,
         -`education`,
         -`political_party`,
         -`region`,
         -`zip`)
race <- race %>% 
  mutate(group = `Respondent.Group`,
         control = case_when(group == 'Control Scenario' ~ 1,  
                             group != 'Control Scenario' ~ 0),
         t1_we = case_when(group == 'Scenario #1, Treatment' ~ 1,
                        group != 'Scenario #1, Treatment' ~ 0),
         t2_wk = case_when(group == 'Scenario #2, Treatment' ~ 1,
                        group != 'Scenario #2, Treatment' ~ 0),
         t3_be = case_when(group == 'Scenario #3, Treatment' ~ 1,
                        group != 'Scenario #3, Treatment' ~ 0),
         t4_bk = case_when(group == 'Scenario #4, Treatment' ~ 1,
                        group != 'Scenario #4, Treatment' ~ 0),
         t5_ire = case_when(group == 'Scenario #5, Treatment' ~ 1,
                        group != 'Scenario #5, Treatment' ~ 0),
         t6_irk = case_when(group == 'Scenario #6, Treatment' ~ 1,
                        group != 'Scenario #6, Treatment' ~ 0))
race <- race %>% 
  select(-`Respondent.Group`) 
race <- race %>% 
  mutate(skin = ifelse(control == 1, "Control",
                       ifelse(t1_we == 1 | t2_wk == 1, 'White', 
                       ifelse(t3_be == 1 | t4_bk == 1, 'Black', 'IR'))))
race <- race %>% 
  mutate(geography = ifelse(control == 1, "Control",
                            ifelse(t1_we == 1 | t3_be == 1 | t5_ire == 1, 'Estonia', 'Kuwait'))) 

race <- race %>% 
  mutate(skin_binary = ifelse(control == 1, "Control",
                               ifelse(t5_ire == 1 | t6_irk == 1, "NA",
                              ifelse(t1_we == 1 | t2_wk == 1, 'White', 'Dark')))) 
                                    
race <- race %>% 
  mutate(geography_binary = ifelse(control == 1, "Control",
                                   ifelse(t1_we == 1 | t3_be == 1 |  t5_ire == 1, 'Western', 'Non-Western')))

race <- race %>% 
  rename(time = Duration..in.seconds., sex = Q1, age = Q2, ethnicity = Q3, ed = Q4, check = Q5, supt = Q6, supt_time = Q6A_Page.Submit, leg = Q7, leg_time = Q7A_Page.Submit, responsible = Q8_1, responsible_time = Q8A_Page.Submit, manipulation1 = Q9, manipulation2 = Q10, open_ended = Q11, t1_party = Q12, t1_ethno1 = Q13_1, t1_ethno2 = Q13_2, t1_ethno3 = Q13_3, t1_ethno4 = Q13_4, t1_ethno5 = Q13_5, t1_ethno6 = Q13_6, t1_ethno7 = Q13_7, t1_force = Q14, t1_internationalism = Q15, t1_religion = Q16, t1_religiosity = Q17, t1_morality = Q18, t1_congress = Q19, t1_legality = Q20, t1_emotion1 = Q21_1, t1_emotion2 = Q21_2, t1_emotion3 = Q21_3, t1_emotion4 = Q21_4, t1_emotion5 = Q21_5, t1_emotion6 = Q21_6, t1_emotion7 = Q21_7, t1_emotion8 = Q21_8, t1_emotion9 = Q21_9, t1_race1 = Q22_1, t1_race2 = Q22_2, t1_race3 = Q22_3, t1_race4 = Q22_4, t1_race5 = Q22_5, t1_race6 = Q22_6)
race <- race %>% 
  rename(t2_party = Q23, t2_ethno1 = Q24_1, t2_ethno2 = Q24_2, t2_ethno3 = Q24_3, t2_ethno4 = Q24_4, t2_ethno5 = Q24_5, t2_ethno6 = Q24_6, t2_ethno7 = Q24_7, t2_force = Q25, t2_internationalism = Q26, t2_religion = Q27, t2_religiosity = Q28, t2_morality = Q29, t2_congress = Q30, t2_legality = Q31, t2_emotion1 = Q32_1, t2_emotion2 = Q32_2, t2_emotion3 = Q32_3, t2_emotion4 = Q32_4, t2_emotion5 = Q32_5, t2_emotion6 = Q32_6, t2_emotion7 = Q32_7, t2_emotion8 = Q32_8, t2_emotion9 = Q32_9, t2_race1 = Q33_1, t2_race2 = Q33_2, t2_race3 = Q33_3, t2_race4 = Q33_4, t2_race5 = Q33_5, t2_race6 = Q33_6)
race <- race %>% 
  rename(t3_party = Q34, t3_ethno1 = Q35_1, t3_ethno2 = Q35_2, t3_ethno3 = Q35_3, t3_ethno4 = Q35_4, t3_ethno5 = Q35_5, t3_ethno6 = Q35_6, t3_ethno7 = Q35_7, t3_force = Q36, t3_internationalism = Q37, t3_religion = Q38, t3_religiosity= Q39, t3_morality = Q40, t3_congress = Q41, t3_legality = Q42, t3_emotion1 = Q43_1, t3_emotion2 = Q43_2, t3_emotion3 = Q43_3, t3_emotion4 = Q43_4, t3_emotion5 = Q43_5, t3_emotion6 = Q43_6, t3_emotion7 = Q43_7, t3_emotion8 = Q43_8, t3_emotion9 = Q43_9, t3_race1 = Q44_1, t3_race2 = Q44_2, t3_race3 = Q44_3, t3_race4 = Q44_4, t3_race5 = Q44_5, t3_race6 = Q44_6)
race <- race %>% 
  rename(t4_party = Q45, t4_ethno1 = Q46_1, t4_ethno2 = Q46_2, t4_ethno3 = Q46_3, t4_ethno4 = Q46_4, t4_ethno5 = Q46_5, t4_ethno6 = Q46_6, t4_ethno7 = Q46_7, t4_force = Q47, t4_internationalism = Q48, t4_religion = Q49, t4_religiosity= Q50, t4_morality = Q51, t4_congress = Q52, t4_legality = Q53, t4_emotion1 = Q54_1, t4_emotion2 = Q54_2, t4_emotion3 = Q54_3, t4_emotion4 = Q54_4, t4_emotion5 = Q54_5, t4_emotion6 = Q54_6, t4_emotion7 = Q54_7, t4_emotion8 = Q54_8, t4_emotion9 = Q54_9, t4_race1 = Q55_1, t4_race2 = Q55_2, t4_race3 = Q55_3, t4_race4 = Q55_4, t4_race5 = Q55_5, t4_race6 = Q55_6)
race <- race %>% 
  rename(t5_party = Q56, t5_ethno1 = Q57_1, t5_ethno2 = Q57_2, t5_ethno3 = Q57_3, t5_ethno4 = Q57_4, t5_ethno5 = Q57_5, t5_ethno6 = Q57_6, t5_ethno7 = Q57_7, t5_force = Q58, t5_internationalism = Q59, t5_religion = Q60, t5_religiosity= Q61, t5_morality = Q62, t5_congress = Q63, t5_legality = Q64, t5_emotion1 = Q65_1, t5_emotion2 = Q65_2, t5_emotion3 = Q65_3, t5_emotion4 = Q65_4, t5_emotion5 = Q65_5, t5_emotion6 = Q65_6, t5_emotion7 = Q65_7, t5_emotion8 = Q65_8, t5_emotion9 = Q65_9, t5_race1 = Q66_1, t5_race2 = Q66_2, t5_race3 = Q66_3, t5_race4 = Q66_4, t5_race5 = Q66_5, t5_race6 = Q66_6)
race <- race %>% 
  rename(t6_party = Q67, t6_ethno1 = Q68_1, t6_ethno2 = Q68_2, t6_ethno3 = Q68_3, t6_ethno4 = Q68_4, t6_ethno5 = Q68_5, t6_ethno6 = Q68_6, t6_ethno7 = Q68_7, t6_force = Q69, t6_internationalism = Q70, t6_religion = Q71, t6_religiosity= Q72, t6_morality = Q73, t6_congress = Q74, t6_legality = Q75, t6_emotion1 = Q76_1, t6_emotion2 = Q76_2, t6_emotion3 = Q76_3, t6_emotion4 = Q76_4, t6_emotion5 = Q76_5, t6_emotion6 = Q76_6, t6_emotion7 = Q76_7, t6_emotion8 = Q76_8, t6_emotion9 = Q76_9, t6_race1 = Q77_1, t6_race2 = Q77_2, t6_race3 = Q77_3, t6_race4 = Q77_4, t6_race5 = Q77_5, t6_race6 = Q77_6) 
race <- race %>% 
  rename(ctl_party = Q78, ctl_ethno1 = Q79_1, ctl_ethno2 = Q79_2, ctl_ethno3 = Q79_3, ctl_ethno4 = Q79_4, ctl_ethno5 = Q79_5, ctl_ethno6 = Q79_6, ctl_ethno7 = Q79_7, ctl_force = Q80, ctl_internationalism = Q81, ctl_religion = Q82, ctl_religiosity  = Q83, ctl_morality = Q84, ctl_congress = Q85, ctl_legality = Q86, ctl_emotion1 = Q87_1, ctl_emotion2 = Q87_2, ctl_emotion3 = Q87_3, ctl_emotion4 = Q87_4, ctl_emotion5 = Q87_5, ctl_emotion6 = Q87_6, ctl_emotion7 = Q87_7, ctl_emotion8 = Q87_8, ctl_emotion9 = Q87_9, ctl_race1 = Q88_1, ctl_race2 = Q88_2, ctl_race3 = Q88_3, ctl_race4 = Q88_4, ctl_race5 = Q88_5, ctl_race6 = Q88_6) 

race <- race %>% 
  rename(racist = Q89, income = Q90, state = Q91, military = Q92, years = Q93, branch = Q94, sof = Q95, rank = Q96, deployments = Q97, combat = Q98)
race <- race %>% mutate(across(!open_ended & !group & !skin & !geography & !skin_binary & !geography_binary, as.numeric))
```

### 1) Summary Statistics: 

``````{r summary, echo=FALSE, eval=TRUE, message=FALSE}
race_table <- race %>% select(sex, age, ethnicity, ed, income, group) %>% filter(sex != 3) %>% filter(group != 'Control Scenario')
race_table <- race_table %>% rename(Sex = sex, Age = age, Ethnicity = ethnicity, Education = ed, Income = income)
race_table$Sex <- ifelse(race_table$Sex==1, "Men", 
                         ifelse(race_table$Sex==2, "Women", NA))
race_table$Age <- ifelse(race_table$Age==1, "Under 18", 
                          ifelse(race_table$Age==2, "19-25",
                                 ifelse(race_table$Age==3, "26-35", 
                                        ifelse(race_table$Age==4, "36-45",
                                               ifelse(race_table$Age==5, "46-55", 
                                                      ifelse(race_table$Age==6, "56-65",
                                                             ifelse(race_table$Age==7, "Over 66",NA)))))))

race_table$Ethnicity <- ifelse(race_table$Ethnicity==1, "American Indian, Alaskan Native", 
                          ifelse(race_table$Ethnicity==2, "Asian",
                                 ifelse(race_table$Ethnicity==3, "Black", 
                                        ifelse(race_table$Ethnicity==4, "Hispanic",
                                               ifelse(race_table$Ethnicity==5, "Native Hawaiian, Other Pacific Islander", 
                                                      ifelse(race_table$Ethnicity==6, "White, Non-Hispanic", 
                                                             ifelse(race_table$Ethnicity==7, "Other",NA)))))))
race_table$Education <- ifelse(race_table$Education==1, "Some High School", 
                          ifelse(race_table$Education==2, "High School",
                                 ifelse(race_table$Education==3, "Some College", 
                                        ifelse(race_table$Education==4, "2-Year Degree",
                                               ifelse(race_table$Education==5, "4-Year Degree", 
                                                      ifelse(race_table$Education==6, "Advanced Degree", NA))))))
race_table$Income <- ifelse(race_table$Income==1, "< $9,999", 
                          ifelse(race_table$Income==2, "$10,000-$24,999",
                                 ifelse(race_table$Income==3, "$25,000-$49,999", 
                                        ifelse(race_table$Income==4, "$50,000-$74,999",
                                               ifelse(race_table$Income==5, "$75,000-$99,999", 
                                                      ifelse(race_table$Income==6, "> $100,000", NA))))))
race_table$group <- ifelse(race_table$group=='Control Scenario', "CTL", 
                          ifelse(race_table$group=='Scenario #1, Treatment', "T1 (WH/EE)",
                                 ifelse(race_table$group=='Scenario #2, Treatment', "T2 (WH/KU)",
                                        ifelse(race_table$group=='Scenario #3, Treatment', "T3 (BL/EE)",
                                               ifelse(race_table$group=='Scenario #4, Treatment', "T4 (BL/KU)",
                                                      ifelse(race_table$group=='Scenario #5, Treatment', "T5 (IR/EE)",
                                                              ifelse(race_table$group=='Scenario #6, Treatment', "T6 (IR/KU)", NA )))))))

theme_gtsummary_journal(journal = "qjecon")
theme_gtsummary_compact()


race_table %>%
  mutate(Sex = factor(Sex, levels = c ("Men", 
                                       "Women")))%>%
  mutate(Age = factor(Age, levels = c ("Under 18", 
                                       "19-25",
                                       "26-35",
                                       "36-45",
                                       "46-55",
                                       "56-65",
                                       "Over 66"))) %>%
  mutate(Ethnicity = factor(Ethnicity, levels = c ("American Indian, Alaskan Native", 
                                                  "Asian",
                                                  "Black",
                                                  "Hispanic", 
                                                  "Native Hawaiian, Other Pacific Islander", 
                                                  "White, Non-Hispanic",
                                                  "Other"))) %>% 
  mutate(Education = factor(Education, levels = c ("Some High School", 
                                                  "High School",
                                                  "Some College",
                                                  "2-Year Degree", 
                                                  "4-Year Degree", 
                                                  "Advanced Degree"))) %>% 
  mutate(Income = factor(Income, levels = c("< $9,999", 
                                            "$10,000-$24,999",
                                            "$25,000-$49,999",
                                            "$50,000-$74,999", 
                                            "$75,000-$99,999", 
                                            "> $100,000"))) %>% 
  tbl_summary(by = group) %>% 
  add_overall() %>%
  modify_header(label ~ "**Group**") %>%
  modify_footnote(update = everything() ~ NA) %>%
  bold_labels
```

### 2) Public Support by Experimental Groups: 

#### a) T-Tests for Experimental Groups: 

```{r support ttests1, echo=FALSE, eval=TRUE }
ctl <- race %>% filter(control == 1) 
t1 <- race %>% filter(t1_we == 1) 
t2 <- race %>% filter(t2_wk == 1) 
t3 <- race %>% filter(t3_be == 1) 
t4 <- race %>% filter(t4_bk == 1)
t5 <- race %>% filter(t5_ire == 1)
t6 <- race %>% filter(t6_irk == 1)
t.test(t6$supt,t1$supt)
t.test(t6$supt,t2$supt)
t.test(t6$supt,t3$supt)
t.test(t6$supt,t4$supt)
t.test(t2$supt,t4$supt)
t.test(t5$supt,t1$supt)
t.test(t5$supt,t2$supt)
t.test(t5$supt,t3$supt)
t.test(t6$supt,t4$supt)
t.test(t1$supt,t3$supt)
t.test(t6$supt,t5$supt)
t.test(ctl$supt,t1$supt)
t.test(ctl$supt,t2$supt)
t.test(ctl$supt,t3$supt)
t.test(ctl$supt,t4$supt)
t.test(ctl$supt,t5$supt)
t.test(ctl$supt,t6$supt)
```

#### b) T-Tests between Skin Colors: 

```{r support ttests2, echo=FALSE, eval=TRUE }
black <- race %>% filter(skin_binary == "Dark") 
white <- race %>% filter(skin_binary == "White") 
t.test(black$supt,white$supt)
```

#### c) T-Tests between Countries: 

```{r racism ttests2, echo=FALSE, eval=TRUE }
ee <- race %>% filter(geography == "Estonia") 
ku <- race %>% filter(geography == "Kuwait") 
t.test(ee$supt,ku$supt)
```

#### d) Support Graphs by Country: 

```{r support by treatment groupsKU, echo=FALSE}
support_summaryKU <- support_summary %>%
  mutate(group = case_when(group == "Scenario #1, Treatment" ~ "t1_we",
                           group == "Scenario #2, Treatment" ~ "t2_wk",
                           group == "Scenario #3, Treatment" ~ "t3_be",
                           group == "Scenario #4, Treatment" ~ "t4_bk",
                           group == "Scenario #5, Treatment" ~ "t5_ire",
                           group == "Scenario #6, Treatment" ~ "t6_irk"))
support_summaryKU <- support_summaryKU[-c(1,3,5),]
ggplot(support_summaryKU, aes(group, supt)) +
  #geom_col(fill = "lightgray") +
  geom_pointrange(aes(ymin = supt-ic, ymax = supt+ic), width = 0.2) +
  # geom_errorbar(aes(ymin = supt, ymax = supt+sd), width = 0.2) +
  scale_x_discrete(labels=c("White*",
                            "Black",
                            "IR")) +
  scale_y_continuous(limits = c(3.25,4.25)) +  
  geom_hline(yintercept = mean(t6$supt), linetype = 2) +
  #coord_cartesian(ylim=c(1,5)) + https://stackoverflow.com/questions/32505298/explain-ggplot2-warning-removed-k-rows-containing-missing-values
  theme_bw() +
  ylab("Mean Support for the U.S. Drone Strike-Kuwait") +
  xlab("Experimental Condition") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ggtitle("")
```

```{r support by treatment groupsES, echo=FALSE}
support_summaryES <- support_summary %>%
  mutate(group = case_when(group == "Scenario #1, Treatment" ~ "t1_we",
                           group == "Scenario #2, Treatment" ~ "t2_wk",
                           group == "Scenario #3, Treatment" ~ "t3_be",
                           group == "Scenario #4, Treatment" ~ "t4_bk",
                           group == "Scenario #5, Treatment" ~ "t5_ire",
                           group == "Scenario #6, Treatment" ~ "t6_irk"))
support_summaryES <- support_summaryES[-c(2,4,6),]
ggplot(support_summaryES, aes(group, supt)) +
  #geom_col(fill = "lightgray") +
  geom_pointrange(aes(ymin = supt-ic, ymax = supt+ic), width = 0.2) +
  # geom_errorbar(aes(ymin = supt, ymax = supt+sd), width = 0.2) +
  scale_x_discrete(labels=c("White",
                            "Black",
                            "IR")) +
  scale_y_continuous(limits = c(3.25,4.25)) +  
  geom_hline(yintercept = mean(t5$supt), linetype = 2) +
  #coord_cartesian(ylim=c(1,5)) + https://stackoverflow.com/questions/32505298/explain-ggplot2-warning-removed-k-rows-containing-missing-values
  theme_bw() +
  ylab("Mean Support for the U.S. Drone Strike-Estonia") +
  xlab("Experimental Condition") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ggtitle("")
```

#### e) Support Graph by Geographic Setting: 

```{r isolated support, echo=FALSE, eval=TRUE, include=TRUE}
support_isolated_summary <- race %>% filter(skin != 'Control') %>% 
  group_by(skin, geography) %>%
  summarise(n=n(),
            sd = sd(supt, na.rm = TRUE),
            se = se(supt, na.rm = TRUE),
            ic = se*qt((1-0.05)/2 + .5, n-1),
            supt = mean(supt))
support_isolated_summary
ggplot(support_isolated_summary, aes(geography, supt)) +
  geom_col(aes(fill = skin), position = position_dodge(0.8), width = 0.7) +
  geom_errorbar(aes(ymin = supt-ic, ymax = supt+ic, group = skin), width = 0.2, position = position_dodge(0.8)) +
  scale_fill_manual(values = c("grey90", "grey70", "grey50")) +
  theme_bw() +
  labs(fill = "Skin Color") +
  ylab("Mean Support of the U.S. Drone Strike") +
  xlab("Geographic Setting") +
  # theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ggtitle("")
```
 
### 3) Dataset for Regression Analysis: 

```{r update dataset, echo=FALSE}
race <- race %>%  
  mutate(party_score = ifelse(control == 1, ctl_party,  
                              ifelse(t1_we == 1, t1_party,  
                                     ifelse(t2_wk == 1, t2_party,
                                            ifelse(t3_be == 1, t3_party,
                                                   ifelse(t4_bk == 1, t4_party,
                                                          ifelse(t5_ire == 1, t5_party, 
                                                                 ifelse(t6_irk == 1, t6_party,NA))))))))
race <- race %>%  
  mutate(force_score = ifelse(control == 1, ctl_force,  
                              ifelse(t1_we == 1, t1_force,  
                                     ifelse(t2_wk == 1, t2_force,
                                            ifelse(t3_be == 1, t3_force,
                                                   ifelse(t4_bk == 1, t4_force,
                                                          ifelse(t5_ire == 1, t5_force, 
                                                                 ifelse(t6_irk == 1, t6_force,NA))))))))
race <- race %>%  
  mutate(internationalism_score = ifelse(control == 1, ctl_internationalism,  
                              ifelse(t1_we == 1, t1_internationalism,  
                                     ifelse(t2_wk == 1, t2_internationalism,
                                            ifelse(t3_be == 1, t3_internationalism,
                                                   ifelse(t4_bk == 1, t4_internationalism,
                                                          ifelse(t5_ire == 1, t5_internationalism, 
                                                                 ifelse(t6_irk== 1, t6_internationalism,NA))))))))
race <- race %>%  
  mutate(religiosity_score = ifelse(control == 1, ctl_religiosity ,  
                              ifelse(t1_we == 1, t1_religiosity ,  
                                     ifelse(t2_wk == 1, t2_religiosity , 
                                            ifelse(t3_be == 1, t3_religiosity ,
                                                   ifelse(t4_bk == 1, t4_religiosity ,
                                                          ifelse(t5_ire == 1, t5_religiosity , 
                                                                 ifelse(t6_irk == 1, t6_religiosity ,NA))))))))
race <- race %>%  
  mutate(morality_score = ifelse(control == 1, ctl_morality,  
                              ifelse(t1_we == 1, t1_morality,  
                                     ifelse(t2_wk == 1, t2_morality, 
                                            ifelse(t3_be == 1, t3_morality,
                                                   ifelse(t4_bk == 1, t4_morality,
                                                          ifelse(t5_ire == 1, t5_morality, 
                                                                 ifelse(t6_irk == 1, t6_morality,NA))))))))
race <- race %>%  
  mutate(congress_score = ifelse(control == 1, ctl_congress,  
                              ifelse(t1_we == 1, t1_congress,  
                                     ifelse(t2_wk == 1, t2_congress, 
                                            ifelse(t3_be == 1, t3_congress,
                                                   ifelse(t4_bk == 1, t4_congress,
                                                          ifelse(t5_ire == 1, t5_congress, 
                                                                 ifelse(t6_irk == 1, t6_congress,NA))))))))
race <- race %>%  
  mutate(legality_score = ifelse(control == 1, ctl_legality,  
                              ifelse(t1_we == 1, t1_legality,  
                                     ifelse(t2_wk == 1, t2_legality, 
                                            ifelse(t3_be == 1, t3_legality,
                                                   ifelse(t4_bk == 1, t4_legality,
                                                          ifelse(t5_ire == 1, t5_legality, 
                                                                 ifelse(t6_irk == 1, t6_legality,NA))))))))
race <- race %>% 
  mutate(ethno_score = ifelse(control == 1, ((ctl_ethno1+ctl_ethno2+ctl_ethno3+ctl_ethno4+ctl_ethno5+ctl_ethno6+ctl_ethno7)/7), 
                       ifelse(t1_we == 1, ((t1_ethno1+t1_ethno2+t1_ethno3+t1_ethno4+t1_ethno5+t1_ethno6+t1_ethno7)/7), 
                       ifelse(t2_wk == 1, ((t2_ethno1+t2_ethno2+t2_ethno3+t2_ethno4+t2_ethno5+t2_ethno6+t2_ethno7)/7), 
                       ifelse(t3_be == 1, ((t3_ethno1+t3_ethno2+t3_ethno3+t3_ethno4+t3_ethno5+t3_ethno6+t3_ethno7)/7), 
                       ifelse(t4_bk == 1, ((t4_ethno1+t4_ethno2+t4_ethno3+t4_ethno4+t4_ethno5+t4_ethno6+t4_ethno7)/7),
                       ifelse(t5_ire == 1, ((t5_ethno1+t5_ethno2+t5_ethno3+t5_ethno4+t5_ethno5+t5_ethno6+t5_ethno7)/7), 
                       ifelse(t6_irk == 1, ((t6_ethno1+t6_ethno2+t6_ethno3+t6_ethno4+t6_ethno5+t6_ethno6+t6_ethno7)/7),NA))))))))
race <- race %>% 
  mutate(emotion_score = 
           ifelse(control == 1, ((ctl_emotion1+ctl_emotion2+ctl_emotion3+ctl_emotion4+ctl_emotion5+ctl_emotion6+ctl_emotion7+ctl_emotion8+ctl_emotion9)/9), 
           ifelse(t1_we == 1, ((t1_emotion1+t1_emotion2+t1_emotion3+t1_emotion4+t1_emotion5+t1_emotion6+t1_emotion7+t1_emotion8+t1_emotion9)/9), 
           ifelse(t2_wk == 1, ((t2_emotion1+t2_emotion2+t2_emotion3+t2_emotion4+t2_emotion5+t2_emotion6+t2_emotion7+t2_emotion8+t2_emotion9)/9),
           ifelse(t3_be == 1, ((t3_emotion1+t3_emotion2+t3_emotion3+t3_emotion4+t3_emotion5+t3_emotion6+t3_emotion7+t3_emotion8+t3_emotion9)/9),
           ifelse(t4_bk == 1, ((t4_emotion1+t4_emotion2+t4_emotion3+t4_emotion4+t4_emotion5+t4_emotion6+t4_emotion7+t4_emotion8+t4_emotion9)/9),
           ifelse(t5_ire == 1, ((t5_emotion1+t5_emotion2+t5_emotion3+t5_emotion4+t5_emotion5+t5_emotion6+t5_emotion7+t5_emotion8+t5_emotion9)/9), 
           ifelse(t6_irk == 1, ((t6_emotion1+t6_emotion2+t6_emotion3+t6_emotion4+t6_emotion5+t6_emotion6+t6_emotion7+t6_emotion8+t6_emotion9)/9),NA))))))))
race <- race %>% 
  mutate(race_score = 
           ifelse(control == 1, ((ctl_race1+ctl_race2+ctl_race3+ctl_race4+ctl_race5+ctl_race6)), 
           ifelse(t1_we == 1, ((t1_race1+t1_race2+t1_race3+t1_race4+t1_race5+t1_race6)), 
           ifelse(t2_wk == 1, ((t2_race1+t2_race2+t2_race3+t2_race4+t2_race5+t2_race6)),
           ifelse(t3_be == 1, ((t3_race1+t3_race2+t3_race3+t3_race4+t3_race5+t3_race6)),
           ifelse(t4_bk == 1, ((t4_race1+t4_race2+t4_race3+t4_race4+t4_race5+t4_race6)),
           ifelse(t5_ire == 1, ((t5_race1+t5_race2+t5_race3+t5_race4+t5_race5+t5_race6)), 
           ifelse(t6_irk == 1, ((t6_race1+t6_race2+t6_race3+t6_race4+t6_race5+t6_race6)),NA))))))))
race_regress <- race %>% select(time,sex,age,ed,supt,ethnicity,income,military,group,ethno_score,emotion_score,party_score,force_score,
                                internationalism_score,religiosity_score,morality_score,congress_score,legality_score,racist, race_score, skin,
                                geography,skin_binary, geography_binary, manipulation1, manipulation2,control,t1_we,t2_wk,t3_be,t4_bk,t5_ire,
                                t6_irk) %>% 
                                filter(party_score != 8)
```

### 4) Regression Analysis:

#### a) Kuwait: 

```{r results='asis', echo=FALSE, eval=TRUE}
race_regressKU <- race_regress
race_regressKU %<>% mutate(new_group = ifelse(t1_we == 1, 'Yes', 
                                             ifelse(t3_be == 1, 'Yes',
                                                     ifelse(t5_ire == 1, 'Yes', 'No')))) %<>% filter(new_group == 'No')
race_regressKU %<>% 
  mutate(group=factor(group)) %<>% 
  mutate(group=fct_relevel(group,c("Scenario #6, Treatment","Scenario #2, Treatment", "Scenario #4, Treatment")))
race_regressKU2 <- race_regressKU 
race_regressKU2 %<>% mutate(overall = (time)/60) %>% 
  filter(overall < 6.682614 & overall > 1)
race_regressKU3 <- race_regressKU
race_regressKU3 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
race_regressKU4 <- race_regressKU2
race_regressKU4 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
model1 <- lm(supt~group,data=race_regressKU)
model2 <- lm(supt~group+sex+age+income+ethnicity,data=race_regressKU)
model3 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressKU)
model4 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressKU2)
model5 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressKU3)
model6 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressKU4)
stargazer(model1,model2,model3,model4,model5,model6,type = "html", style = "qje",
          digits = 2,
          ci = TRUE,
          title = "Public Support for U.S. Counterterrorism Drone Strikes", 
          covariate.labels = c("White, Kuwait","Black, Kuwait", "Sex", "Age", "Income", "Race", "Ideology","Ethnocentrism","Use of Force","Morality","Domestic 
          Approval","International Law","Racial Resentment","Strike are Racist","Military Service"), 
          dep.var.labels   = "Support",
          df = FALSE,
          omit.stat = c("rsq", "ser")) #regression models for Appendix C, SI
```

#### b) Estonia: 

```{r results='asis', echo=FALSE, eval=TRUE}
race_regressES <- race_regress
race_regressES %<>% mutate(new_group = ifelse(t1_we == 1, 'Yes', 
                                             ifelse(t3_be == 1, 'Yes',
                                                     ifelse(t5_ire == 1, 'Yes', 'No')))) %<>% filter(new_group == 'Yes')
race_regressES %<>% 
  mutate(group=factor(group)) %<>% 
  mutate(group=fct_relevel(group,c("Scenario #5, Treatment","Scenario #1, Treatment", "Scenario #3, Treatment")))
race_regressES2 <- race_regressES
race_regressES2 %<>% mutate(overall = (time)/60) %>% 
  filter(overall < 6.682614 & overall > 1)
race_regressES3 <- race_regressES
race_regressES3 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
race_regressES4 <- race_regressES2
race_regressES4 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
model7 <- lm(supt~group,data=race_regressES)
model8 <- lm(supt~group+sex+age+income+ethnicity,data=race_regressES)
model9 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressES)
model10 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressES2)
model11 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressES3)
model12 <- lm(supt~group+sex+age+income+ethnicity+party_score+ethno_score+force_score+morality_score+congress_score+legality_score+race_score+racist+military,data=race_regressES4)
stargazer(model7,model8,model9,model10,model11,model12,type = "html", style = "qje",
          digits = 2,
          ci = TRUE,
          title = "Public Support for U.S. Counterterrorism Drone Strikes", 
          covariate.labels = c("White, Estonia","Black, Estonia", "Sex", "Age", "Income", "Race", "Ideology","Ethnocentrism","Use of Force","Morality","Domestic 
          Approval","International Law","Racial Resentment","Strike are Racist","Military Service"), 
          dep.var.labels   = "Support",
          df = FALSE,
          omit.stat = c("rsq", "ser")) #regression models for Appendix C, SI 
```


### 5) Robustness Check for Manipulation Check 

#### a) Kuwait:

```{r robustness check for manipulation checkKU, echo=FALSE}
support_summary2 <- race %>% filter (group != "Control Scenario") 
support_summary2 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
t2 <- support_summary2 %>% filter(t2_wk == 1) 
t4 <- support_summary2 %>% filter(t4_bk == 1)
t6 <- support_summary2 %>% filter(t6_irk == 1)
t.test(t6$supt,t2$supt)
t.test(t6$supt,t4$supt)
t.test(t2$supt,t4$supt)
support_summary2 %<>% filter (group != "Control Scenario") %>%
  group_by(group) %>%
  summarise(n=n(),
            sd = sd(supt, na.rm = TRUE),
            se = se(supt, na.rm = TRUE),
            ic = se*qt((1-0.05)/2 + .5, n-1),
            supt = mean(supt))
support_summaryKU2 <- support_summary2 %>%
  mutate(group = case_when(group == "Scenario #1, Treatment" ~ "t1_we",
                           group == "Scenario #2, Treatment" ~ "t2_wk",
                           group == "Scenario #3, Treatment" ~ "t3_be",
                           group == "Scenario #4, Treatment" ~ "t4_bk",
                           group == "Scenario #5, Treatment" ~ "t5_ire",
                           group == "Scenario #6, Treatment" ~ "t6_irk"))
support_summaryKU2 <- support_summaryKU2[-c(1,3,5),]
ggplot(support_summaryKU2, aes(group, supt)) +
  #geom_col(fill = "lightgray") +
  geom_pointrange(aes(ymin = supt-ic, ymax = supt+ic), width = 0.2) +
  # geom_errorbar(aes(ymin = supt, ymax = supt+sd), width = 0.2) +
  scale_x_discrete(labels=c("White",
                            "Black**",
                            "IR")) +
  scale_y_continuous(limits = c(3.00,4.25)) +  
  geom_hline(yintercept = mean(t6$supt), linetype = 2) +
  #coord_cartesian(ylim=c(1,5)) + https://stackoverflow.com/questions/32505298/explain-ggplot2-warning-removed-k-rows-containing-missing-values
  theme_bw() +
  ylab("Mean Support for the U.S. Drone Strike-Kuwait") +
  xlab("Experimental Condition") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ggtitle("")
```

#### b) Estonia:

```{r robustness check for manipulation checkKU, echo=FALSE}
support_summary3 <- race %>% filter (group != "Control Scenario") 
support_summary3 %<>% filter(skin != 'Control') %<>%
  mutate(new_man = ifelse(skin == 'Black' & manipulation1 == 1, 'Yes',
                                             ifelse(skin == 'White' & manipulation1 == 2, 'Yes', 
                                                  ifelse(skin == 'IR' & manipulation1 == 3, 'Yes','No')))) %<>% filter(new_man == 'Yes')
t1 <- support_summary3 %>% filter(t1_we == 1) 
t3 <- support_summary3 %>% filter(t3_be == 1)
t5 <- support_summary3 %>% filter(t5_ire == 1)
t.test(t5$supt,t1$supt)
t.test(t5$supt,t3$supt)
t.test(t1$supt,t3$supt)
support_summary3 %<>% filter (group != "Control Scenario") %>%
  group_by(group) %>%
  summarise(n=n(),
            sd = sd(supt, na.rm = TRUE),
            se = se(supt, na.rm = TRUE),
            ic = se*qt((1-0.05)/2 + .5, n-1),
            supt = mean(supt))
support_summaryES2 <- support_summary3 %>%
  mutate(group = case_when(group == "Scenario #1, Treatment" ~ "t1_we",
                           group == "Scenario #2, Treatment" ~ "t2_wk",
                           group == "Scenario #3, Treatment" ~ "t3_be",
                           group == "Scenario #4, Treatment" ~ "t4_bk",
                           group == "Scenario #5, Treatment" ~ "t5_ire",
                           group == "Scenario #6, Treatment" ~ "t6_irk"))
support_summaryES2 <- support_summaryES2[-c(2,4,6),]
ggplot(support_summaryES2, aes(group, supt)) +
  #geom_col(fill = "lightgray") +
  geom_pointrange(aes(ymin = supt-ic, ymax = supt+ic), width = 0.2) +
  # geom_errorbar(aes(ymin = supt, ymax = supt+sd), width = 0.2) +
  scale_x_discrete(labels=c("White",
                            "Black",
                            "IR")) +
  scale_y_continuous(limits = c(3.00,4.25)) +  
  geom_hline(yintercept = mean(t5$supt), linetype = 2) +
  #coord_cartesian(ylim=c(1,5)) + https://stackoverflow.com/questions/32505298/explain-ggplot2-warning-removed-k-rows-containing-missing-values
  theme_bw() +
  ylab("Mean Support for the U.S. Drone Strike-Estonia") +
  xlab("Experimental Condition") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  ggtitle("")
```